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Abstract. Motivated by the possibility of observing gravitational waves from 
merging black holes whose spins are nearly extremal (i.e., 1 in dimensionless 
units), we present numerical waveforms from simulations of merging black holes 
with the highest spins simulated to date: (1) a 25.5-orbit inspiral, merger, and 
ringdown of two holes with equal masses and spins of magnitude 0.97 aligned with 
the orbital angular momentum; and (2) a previously reported 12.5-orbit inspiral, 
merger, and ringdown of two holes with equal masses and spins of magnitude 
0.95 anti-aligned with the orbital angular momentum. First, we consider the 
horizon mass and spin evolution of the new aligned-spin simulation. During 
the inspiral, the horizon area and spin evolve in remarkably close agreement with 
Alvi's analytic predictions, and the remnant hole's final spin agrees reasonably well 
with several analytic predictions. We also find that the total energy emitted by a 
real astrophysical system with these parameters — almost all of which is radiated 
during the time included in this simulation — would be 10.952% of the initial mass 
at infinite separation. Second, we consider the gravitational waveforms for both 
simulations. After estimating their uncertainties, we compare the waveforms 
to several post-Newtonian approximants, finding significant disagreement well 
before merger, although the phase of the TaylorT4 approximant happens to 
agree remarkably well with the numerical prediction in the aligned-spin case. 
We find that the post-Newtonian waveforms have sufficient uncertainty that 
hybridized waveforms will require far longer numerical simulations (in the absence 
of improved post-Newtonian waveforms) for accurate parameter estimation of low- 
mass binary systems. 



PACS numbers: 04.25.dg, 04.30.-w 



1. Introduction 

In the next decade, advanced ground-based detectors such as the advanced Laser 
Interferometer Gravitational- Wave Observatory (advanced LIGO) [1,2], Virgo [3], and 
the Large-scale Cryogenic Gravitational-wave Telescope (LCGT) [4] are expected to 
directly observe gravitational waves for the first time; coalescing black holes are among 
the most important sources of gravitational waves for these detectors. Numerical 
predictions of binary-black-hole (BBH) waveforms are crucial tools for detecting 
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these waves: for example, the Numerical INJection Analysis (NINJA) project [5, 6] 
is testing gravitational-wave search pipelines using numerical BBH waveforms, and 
the Numerical-Relativity and Analytical-Relativity (NR-AR) project [7] is working 
to calibrate analytic template banks for gravitational-wave searches using numerical 
BBH waveforms. Numerical BBH waveforms are also important tools for parameter 
estimation [8, 9, 10, 11]. 

Beginning with Pretorius's 2005 breakthrough [12], several groups have 
successfully completed numerical simulations of the inspiral, merger, and ringdown 
of BBHs in a variety of configurations (see references [13, 14] for recent reviews). 
Simulations of BBHs with merging holes whose spins are nearly extremal (i.e., ~ 1 in 
dimensionless units, the theoretical upper limit for stationary holes) are a challenging 
but potentially important case, since black holes with nearly extremal spin might 
exist [15, 16, 17, 18, 19, 20, 21, 22] and thus might be among the BBHs emitting 
gravitational waves. Almost all published BBH simulations to date start with initial 
data in which 3 of the 4 Einstein constraint equations are solved analytically using the 
solutions of Bowen and York [23, 24]; this choice of initial data limits the black-hole 
dimensionless spins to x ^ 0.93 [25, 26, 27]. Dain, Lousto, and Zlochower have closely 
approached this "Bowen- York limit" by evolving an equal-mass BBH with equal spins 
of magnitude x — 0.924 aligned with the orbital angular momentum [28]. Note that 
the Bowen- York limit is actually far from extremal in terms of the physical effects of 
the spin: for example, a black hole with spin x = 0-93 has only 59% of the rotational 
energy of an extremal hole of the same mass. 

By using an alternative method to construct BBH initial data, one can surpass 
the Bowen- York limit. In reference [29], three of us (Lovelace, Scheel, and Szilagyi) 
constructed and evolved (through 12.5 orbits, merger, and ringdown) BBH initial 
data (based on a weighted superposition of two boosted, spinning Kerr-Schild black 
holes [30]) with equal masses and equal spins of magnitude x = 0.94905 anti-aligned 
with the orbital angular momentum. In this paper, we present a new BBH simulation 
(through 25.5 orbits of inspiral, merger, and ringdown) with spins of magnitude 
X = 0.96950 aligned with the orbital angular momentum. These simulations are 
the first to surpass the Bowen- York limit and contain the most nearly extremal black 
holes yet simulated, with the black holes with spin x = 0.94905 and x = 0.96950 
having 65% and 72% as much rotational energy, respectively, as an extremal hole of 
the same mass. 

In this paper, we consider the gravitational waveforms from these two simulations. 
We begin in Sec. 2 by summarizing the numerical methods we use to construct and 
evolve rapidly-spinning BBH initial data and also the methods used to extract and 
extrapolate the gravitational waveforms. In Sec. 3, we examine the horizon mass and 
spin evolution in the new x = 0.96950 simulation. Then, in Sec. 4.1, we examine 
the emitted gravitational waveforms and their accuracy. In Sec. 4.2, we compare the 
numerical waves to several post-Newtonian (pN) approximants. We conclude with a 
brief discussion of the implications of our results in Sec. 5. 

2. Numerical methods 

2.1. Initial data 

To construct BBH initial data with rapid spins, we use the method of reference [30] 
and the references therein: we use a spectral elliptic solver [32] to solve the extended 
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Table 1. Some properties of the BBH simulations considered in this paper. 
The quantity q is the mass ratio, M is the sum of the holes' Christodoulou 
masses, and hole i (where i = A,B) has dimensionless spin xf along the z axis 
(i.e., in the direction of the orbital angular momentum). Also listed are the 
Arnowitt-Deser-Misner (ADM) mass A^adm and angular momentum J%j^i^ (e-g-i 
Eqs. (25)-(26) of reference [30]) as well as the mass and spin of the final black 
hole. The quantities M, Mi, and xf are measured after the initial relaxation and 
junk radiation emission (at times t = lOOOM and t = 500M for S^'^ and S^h'^ , 
respectively). Note that the mass M and spin xf are time-dependent; during the 
inspiral from infinite separation to the start of simulations S^-^^ and , Alvi's 
formulas [3f] predict that the total mass M and spin xf would change by less 
than one part in 10® and by less than one part in 10*, respectively. The ADM 
quantities Madm and Jadm are evaluated at t = 0, and Mjjnal and Xflnal 
measured at the final time of the simulation. All uncertainties are estimated as 
the difference between the quantity at highest and second-highest resolution. 
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Table 2. The initial angular velocity f2o, radial velocity ro/ro, and coordinate 
separation do and the corresponding estimated orbital eccentricity e. The 
parameters do, ro/ro, and Qo are chosen when constructing the initial data; for 
simplicity, only the first five significant figures are shown. 

conformal thin sandwich equations with quasi-equihbrium boundary conditions [33, 34, 
35, 36, 37, 38]. Wc choose free data based on a superposition of two boosted, spinning 
Kerr-Schild black holes, tuning the freely specifiable parameters with a numerical root- 
finding algorithm based on Broyden's method [39, 40] to obtain the desired masses and 
spins. We reduced the orbital eccentricity using the iterative technique of reference [41] 
which is based on fits of the orbital frequency. 

We summarize the two configurations we consider in Tables 1 and 2. 

2.2. Evolution 

We evolve the initial data summarized in the Sec. 2.1 using the Spectral Einstein 
Code (SpEC) using the methods summarized in Sec. Ill of reference [29], which 
extend the techniques of reference [42] and the references therein to accommodate 
BBHs with spins above the Bowen-York limit. The evolution (but not the properties 
of the resulting gravitational waveform) of configuration S'^^ was first reported in 
reference [29]; we present the evolution of configuration S^:^ for the first time here. 
Full details our our methods will be given in a future paper; here, we merely summarize 
our method, highlighting some of the additional techniques necessary to merge . 
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Figure 1. The trajectory and constraint violation in simulation SlJ.^.^- Left 
panel: The trajectories of the centers of the individual apparent horizons. Also 
shown are the individual and common horizons at the final time at which the 
trajectories are shown. Right panel: For each resolution, the constraint violation 
(specifically, the normalized constraint energy defined in Eq. (71) of reference [43]) 
normalized over the entire computational domain. Also note that the horizontal 
scale changes at t/M = 6000. 



As described in reference [29], we excise the singularities inside the black 
holes from our computational domain, using a time-dependent, adaptively adjusted 
coordinate mapping to keep the excision surfaces inside the horizons. Because we do 
not apply boundary conditions on the excision surface, the evolution is only well posed 
if the excision surface is a pure-outflow surface (i.e, if it has no incoming characteristic 
fields). During the inspiral, we enforce this condition by controlling the size of the 
excision surface such that it tracks the size of the apparent horizon; shortly before 
merger (i.e., during the final ~ 1 orbit of evolution and during the final ^ 3 orbits 
in the evolutions of S^'^) we control the characteristic speeds directly by adjusting 
the velocity of the excision surface. During the final ^ 0.25 orbits before merger of 
^0.95 during the final ~ 3 orbits in S^'^ we employed the spectral adaptive mesh 
refinement summarized in reference [29] . We also note that when evolving both S^'^ 
and S^'^ ^ we smoothly change gauge conditions to the dampcd-harmonic condition 
described in reference [42] at the beginning of the evolution instead of shortly before 
merger. 

Because the holes spend more time in a highly dynamical and distorted state, we 
find that merging requires that our coordinate mapping must track the apparent- 
horizon shapes more accurately (i.e., to a higher spherical-harmonic resolution t). 
We also find that we must carefully fine-tune the characteristic speed control to 
balance two competing requirements: (1) that the excision surface have no incoming 
characteristic fields; and (2) that the excision surface remain inside the apparent 
horizon. Here, we do this fine-tuning manually in order to merge S'^^] in the future, 
we plan to employ a method that handles any tuning automatically. Because the 
remnant hole also has a rapid spin [particularly just after it forms (figure 2)], we 
similarly control the horizon shape and characteristic speeds during the ringdown of 
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(with some manual fine-tuning). 

To measure the characteristic speeds on the excision surface with sufficient 
accuracy near merger, just before merger in S^'^ we adopt a computational domain 
where the individual apparent horizons lie within a thin, high-resolution spherical shell 
(instead of within a set of cylindrical subdomains, as evolution S^'^ employed). 

The right panel of figure 1 shows the numerical convergence of our method for 
the more demanding case . The constraints initially grow, then drop as the 
initial burst of spurious gravitational radiation leaves the computational domain. The 
constraints then remain clearly convergent throughout the inspiral. Shortly before 
enabling spectral adaptive mesh refinement, we found it necessary to increase the 
resolution of the inner spheres in the low-resolution run in order to control adequately 
the characteristic speeds; this appears as a discontinuous drop in the low-resolution 
constraint energy. As the evolution approaches merger, the constraint violation grows 
in spite of the spectral adaptive mesh refinement. During ringdown, the constraints 
rapidly drop as the hole relaxes to its final state. As the radiation leaves the grid, 
the constraints drop sharply in the low and medium resolutions but not in the high 
resolution J. 

Finally, we briefiy note the computational cost of these two runs. Because the 
simulation involves higher spins, a very large orbital hangup effect (requiring 
twice as many orbits to merge from the same initial separation as the simulation) , 
and a large amount of time in a regime where the spacetime is highly dynamical, the 
S^'^ simulation turned out to be much more computationally expensive than the 
^0.95 giniulation. Specifically, the high-resolution simulation S^'^ required w 110, 000 
cpu hours (w 120 days of wallclock time). For comparison, the high-resolution 5**^-^5 
simulation required « 20,000 cpu hours (w 20 days of wallclock time). 

2.3. Waveform extraction 

We extract the gravitational waveform h using the Regge- Wheeler- Zerilli formal- 
ism [44, 45, 46, 47] on concentric spheres with radii from roughly r = 100 M to 400 M; 
we then extrapolate§ the waveforms to infinite radius using the method of Boyle and 
Mroue [49] with polynomials of order = 4. The waveforms arc decomposed in the 
standard way [50] as modes hg^rn of spin s = —2 spherical harmonics. We use these 
modes to define the amplitude and phase of the waveforms in the usual way: 

Ai.mir) = \rhi,m{T)/M\ (/if,™ (r) = unwind {arg[/if,m(r)]} , (1) 

where the factor of r/M removes the radial dependence from Ai,m and the unwind 
function removes discontinuities of 27r in the data caused by branch cuts [51]. We also 

I After a common horizon forms, our numerical simulation stops, interpolates onto a new 
computational domain with only a single excised region just inside the common horizon, and then 
continues using this new domain. Thus the numerical resolution used during the "plunge" (i.e., 
just before merger) is different than during the "ringdown" (i.e., just after merger). The lack 
of convergence at late times in figure 1 follows from differences in the medium and high plunge 
resolution. All 3 plunge resolutions required different fine-tuning in order to merge; we presume that 
these differences are responsible for the behavior visible in figure 1. We have verified that for the 
highest plunge resolution, the constraints converge with ringdown resolution, even at late times. 
§ We extract at finite radii and then extrapolate, since our computational domain only extends out 
to some finite radius. Extrapolation is meant to eliminate near-field effects and gauge dependence; 
to verify that the extrapolated waves do not retain any residual gauge dependence, the extrapolated 
waves could be compared with waveforms obtained using Cauchy Characteristic Extraction (CCE; 
see reference [48] and the references therein for details). In the future, we plan to compare CCE 
waveforms with the extrapolated waveforms presented in this paper. 
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Figure 2. The spin evolution in simulation S9 ?^ 



. Left panel: The dimensionlcss 
spin X as a function of time for one of the individual horizons. The inset zooms 
in on the spin during the inspiral. Right panel: The dimensionless spin x of the 
common horizon as a function of time. Also shown, as individual data points with 
error bars where applicable, are the final spin predictions given by the analytic 
formulas in references [52] ("A"), [53] ("B"), [20] ("C"), [54] ("D"), [55] ("E"), 
and [56] ("F"). The inset zooms in on the final spin of the common horizon. The 
predictions arc shown at different times simply for convenient separation. 



use the frequency uJi^m{T) = dt (t>e.miT). We will occasionally drop the subscripts on 
these quantities, implicitly referring to {£,m) ~ (2,2). 

3. IVIass and spin evolution in the S'^^f'' simulation 

In figure 2, we show the dimensionless spin as a function of time for each resolution of 
the S^f simulat ion. During the initial relaxation, the holes absorb energy, causing 
the dimensionless spin to quickly relax from x = 0.9700 at time t ~ to x = 0.9695 
at time t = lOOOM. 

Similarly to reference [28], we observe a very large orbital hangup [54, 57, 58] 
during the long inspiral: starting from the same initial coordinate separation (and 
thus at approximately the same initial orbital frequency), case S'^^ requires more than 
twice as many orbits to merge as does case S^^^ (compare the left panel of figure 1 
and the top-left panel of figure 3 in reference [29]) and reaches roughly twice the 
orbital frequency. During this long inspiral, the spin remains above x = 0.969 during 
the first 21.5 orbits but then decreases near merger as the spin angular momentum is 
transformed into orbital angular momentum via tidal interactions. 

The mechanism by which this transformation of angular momentum takes place 
has been described by numerous authors [59, 60, 61, 62, 63] including Alvi [31], 
who gave expressions for the rate at which energy and angular momentum would be 
transferred in comparable-mass binaries. In figure 3 we plot those rates as measured 
in the S'^^ simulation and compare to those predicted by Alvi. Alvi's expression uses 
the post-Newtonian velocity parameter v which we set to v = (Alfiy^^, where fl is 
the orbital angular frequency measured in the simulation. Though this comparison 
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Figure 3. Change in irreducible mass and angular momentum of the horizons 
in the system compared with Alvi's predictions [31]. Left panel: rate of 

change of irreducible (apparent-horizon) mass of the individual horizons. Right 
panel: rate of change of spin of the individual horizons. In each case, the noisier 
(thick gray) line represents data measured directly in the simulation, while the 
smoother (thin black) line is given by Alvi's equation (11). The horizontal axis 
is parameterized by t! = {M il)^/^ , where Q is the orbital frequency measured in 
the simulation. Note that v > 0.5 corresponds to the last ~ 40M of time before 
merger. 

is gauge dependent, we find very good agreement — witliin tlie numerical uncertainty 
until very late in the simulation. A similar comparison for the case is not as 

useful because the numerical uncertainties are far larger, though the result is consistent 
within the larger uncertainties. Note that this transfer of angular momentum is a 
2.5-pN spin effect which is incorporated into the calculation of the pN waveforms in 
section 4.2. 

When the common apparent horizon first forms, its spin is very nearly extremal 
but then quickly relaxes as the common horizon expands, eventually settling to a 
final spin Xfinai = 0.94496 ± 0.00001, which is roughly consistent with the predictions 
of several analytic approximations (right panel of figure 2). The mass of the final 
hole is Mfi„ai/M = 0.89048 ± 0.00002. Alvi's formulas suggest that the mass would 
change by less than a part in 10^ prior to the beginning of our simulation if the binary 
had inspiraled from infinite separation; therefore, a "real" binary would have radiated 
E,^i/M = 1 - Mfine^i/M = 10.952% ± 0.002% of its initial mass (= 12.299% ± 0.003% 
of its final mass) throughout its entire inspiral, merger, and ringdown. This efficiency 
is comparable to that of a supernova (~ 15% of the final core mass radiated; see, 
e.g., equation (18.1.1) of reference [64]) but corresponds to a larger total energy 
radiated (since the mass of a BBH is typically larger than the final core mass after 
a supernova). For comparison, the simulation in reference [65] implies that an equal- 
mass nonspinning binary system would radiate about 5% of its mass, while table 1 
shows that the S^'^ system radiates about 3.17% of its mass. When the merger occurs 
at a frequency in the sensitive band of a gravitational-wave detector, we can therefore 
expect that an aligned-spin system should have significantly larger SNR than a similar 
system with anti-aligned spins. 
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4. Gravitational waveforms and post-Newtonian comparisons 

To compare two waveforms, A and we need to align them by fixing the arbitrary 
relative time and phase offsets. Here A and B may refer to two numerical waveforms 
with different resolutions or extrapolation orders, or A and B may refer to pN and 
numerical waveforms. Following reference [66], we align the waveforms by minimizing 
the difference in their phases over a certain range. Specifically, we minimize || the 
quantity 

S(At, A0) = r [d,Ait) ^^B{t- At) - A(j>f dt . (2) 

Each mode of waveform B is then transformed as 

hi^^{t)^heMt + ^t)o-''^^^^^ ■ (3) 

Note that (j) refers to the phase of the {£, m) = (2, 2) mode only; the values of At and 
A(/) are determined once, then each mode is transformed by this equation. 

The optimal values of At and A^i determined by minimizing equation (2) clearly 
depend on the range of integration (ti,t2)- We choose that range based on the 
frequencies of the waveform [10] so that bj{ti) « 0.033 and ^(^2) ~ 0.038.^ This gives 
us a common basis for comparison of the ranges used in the two cases of aligned and 
anti-aligned spins, despite the very different lengths of time over which they inspiral. 

4-1. Waveform accuracy 

We plot the amplitudes of the three dominant modes and the phase of the dominant 
{i,m) = (2,2) mode in the upper panels of figure 4 (for 55^.+^) and figure 5 (for 
S^^). We also estimate the accuracy of the waveforms by measuring convergence 
with respect to increasing numerical resolution in the simulations and with respect 
to increasing order of the polynomial used for extrapolation to infinite radius. The 
relative amplitude convergence and phase convergence are plotted in the lower panels 
of the two figures. 

The overall uncertainty estimate for a given quantity is the sum of the absolute 
values of the resolution convergence and the extrapolation convergence. That is, in 
the notation of the figures, we estimate 

Uncertainty « |(Med.) - (High)| + |(A^ = 4) - (A^ = 3)| . (4) 

For a waveform to be included in the NINJA-2 data set [5], the amplitude and phase of 
the (2, 2) mode must be accurate at merger to within 5% and 0.5 rad, respectively. The 
case exceeds these requirements. The case, however, exceeds the amplitude 
requirement but does not meet the phase accuracy requirement. Through the time 
of merger, the amplitude uncertainty never exceeds 2%, but the phase uncertainty at 
the time of merger is 0.9 rad. 

While the complete waveform does not meet the NINJA-2 accuracy 

requirement for phase, this simulation's 25.5 orbits far exceed the NINJA-2 length 
requirement of 5 orbits aligned with pN before a frequency of Af uj — 0.075. Indeed, if 

II For fixed At, the optimal A<j) can be obtained analytically, reducing the minimization to a one- 
dimensional problem. See reference [66] for details. 

^ This range corresponds to Sui/uj as 15%, which is somewhat larger than the 10% minimum 
recommended by MacDonald et al. [10]. We use a larger range to ensure that the alignment is 
not skewed by small oscillations in the data. 
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Figure 4. The extrapolated gravitational-wave amplitudes and phase for 
simulation S'^^ . Left panel: The dominant wave amplitudes ^2,2, ^3,2, 
and j44^4 at high resolution (top) and relative differences |i5j42, 2/^2,2! between 
resolutions or between extrapolation orders (bottom). Right panel: The phase 
at high resolution (top) and differences \5<j)\ between resolutions (bottom). When 
computing differences, the waveforms are aligned as in (2) between (t — r,)/M = 
1322 and {t - rt)/M = 2852. The merger time (t - r,)/M = 6411 is the time at 
which the (2, 2) amplitude is maximal, denoted by the vertical dashed lines. If the 
waveform were instead truncated to 5 orbits before merger (the NINJA-2 length 
requirement), the amplitude and phase errors would drop significantly (dotted 
lines). Note that the scale on the horizontal axis changes at (t — rt)/M = 6000 
in each plot for improved visibility of the merger and ringdown. 

we omit the first 5000 M of the waveforms, the simulation still easily meets the 
length requirements. Evaluating the errors by aligning between the beginning of the 
shortened waveforms and M uj = 0.075, the convergence measure improves by a factor 
of 10, and the simulation comfortably exceeds the NINJA-2 accuracy requirements 
with a phase uncertainty of 0.1 rad at the time of merger, as shown by the dotted 
lines in the lower panels of figure 4. Thus, we see that using accuracy requirements 
of this sort without regard to the length of the simulation actually creates a perverse 
incentive to produce shorter numerical waveforms, which decreases the accuracy of 
complete "hybridized" waveforms [11]. An exactly analogous situation occurs when 
the criteria require a particular match with respect to a detector noise curve but do 
not stipulate a mass at which that match should be measured; again, the incentive is 
to produce shorter waveforms. 

4-2. Comparison with post- Newtonian approximations 

The post-Newtonian waveform is constructed in two steps: (1) computation of 
the orbital phase of the binary; and (2) computation of the amplitude of the 
waveform using that phase. For nonspinning systems, the formulas needed for those 
computations have been calculated to 3.0-pN order beyond the leading order in 
amplitude and 3.5-pN beyond leading order in phase. Additional spin-orbit and spin- 
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Figure 5. The extrapolated gravitational-wave amplitudes and phase as in 
figure 4 but for the simulation. Here, alignment is performed between 

(t - Tt)/M = 600 and {t - rt)/M = 1640. This simulation is much shorter than 
the aligned case (about half the number of orbits), even though both simulations 
start at the same coordinate separation and consequently at roughly the same 
frequency. Here the anti-aligned spins cause the holes to merge more quickly than 
an analogous non-spinning binary, whereas in case S^'jJ the binary experiences 
a large orbital hangup. Note that the scale on the horizontal axis changes at 
(t — r«)/M = 3000 in each plot for improved visibility of the merger and ringdown. 



spin terms are available to 3.0-pN order in phase and amplitude, + though not all of 
these have yet been expressed in a useful form for generating waveforms. We use 
the expressions given in appendix 1 of [50] for the flux, orbital energy, and tidal 
heating; and the expressions given in equation (9.4) of reference [68] and appendix 
2 of reference [50] for the waveform amplitudes. The sole addition we make is the 
inclusion of a recently published spin-orbit contribution to the flux. In [50], equation 
(A. 13) should be supplemented by adding a term [67] as 

F{v)^F{v) + ^v^%'' (65Jxa + (65 - 68z.)x.,)]}. (5) 

Though formally high in order, this term is large enough in magnitude (for the spins 
we present) to dominate the next-to-leading-order spin contribution to the flux during 
most of the simulations, dominating even the leading-order term several hundred M 
before merger. 

The orbital phase is computed using an energy-balance equation incorporating the 
rate of change of orbital energy and the loss of that energy in the form of tidal heating 
and the gravitational-wave energy flux J-. Various methods exist — referred to as the 
TaylorTl, T2, T3, and T4 approximants [69, 70, 51] — for integrating these equations, 
all of which should be equivalent to the pN order available, in the sense that they differ 
only by higher-order pN terms. This suggests that all four approximants should agree 
with each other and with the numerical waveform, within the uncertainty of the pN 
approximations. We have re-derived each of these approximants using the expressions 
described above and compared the results of each to the numerical waveforms. 



+ See reference [67] and references therein. 




Figure 6. Comparison between numerical and pN data for the system for 

tlie four pN approximants. Left panel: Phase and relative amplitude errors for 
the dominant (£, m) = (2, 2) mode. Right panel: Relative amplitude errors for the 
next two leading modes, {£,m) = (3,2) and (4,4). In all plots, the shaded gray 
region denotes our estimate for the numerical uncertainty, described by equation 4. 
Finely dotted portions of the graphs indicate negative errors; the pN quantity is 
larger than the numerical quantity in these regions. Note that the scale on the 
horizontal axis changes at {t — r*)/!/ = 6000 in each plot for improved visibility 
of the merger and ringdown. 



Figures 6 and 7 show the pN comparisons for and S^^, respectively. Only 

the phase of the (2, 2) mode is shown, because the phases of other modes are integer 
multiples of this to a good degree of accuracy. The gray region in the background of 
each plot shows our uncertainty estimate for the numerical data of the given quantity, 
given by equation (4). and made to be a non-decreasing function of time after the 
beginning of the alignment region. 

One remarkable feature is that the TaylorT4 approximant captures the phase 
surprisingly well for the system (black line in the top left plot of figure 6); it 

agrees with the numerical data within the uncertainty for roughly 3400M after the 
end of the alignment region (in which it is forced to agree with the numerical data) . 
This brings it within 200M of the merger. Contrast that agreement with the other 
approximants, which disagree with the numerical data immediately after (or even 
before) the end of the alignment region. 

Of course, this agreement of TaylorT4 is presumably pure coincidence, as all 
approximants agree with each other within the uncertainty of the pN approximations. 
The same coincidence was found in the equal-mass nonspinning case [51], but has been 
shown not to carry over to systems with other parameters (see, e.g., [58]). Indeed, 
looking at the S*^^ system (figure 7), we see that the TaylorT4 approximant is actually 
the second worst in that case — the sole worse approximant being TaylorTS, which 
has decreasing orbital frequency starting just before the simulation, and is thus an 
especially poor description of the waveform.* We also point out that the system 



* To be more precise, the TaylorTS approximant expresses the orbital frequency of the binary as 
a function of the pN time to coalescence by inverting a power series (which is used directly for 
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(t-r,)/M (t-r,)/M 



Figure 7. Comparison between numerical and pN data as in figure 6 but for tlie 
4^0.95 gystem. Note tfie difference in fiorizontal scale compared to figure 6. Note 
that tfie scale on the horizontal axis changes at (t — rt)/M = 3000 in each plot 
for improved visibility of the merger and ringdowrn. 



inspirals for roughly twice as fong as the S'^^ system and reaches nearly twice the 
frequency, making the phase coherence of TaylorT4 all the more surprising. 

The figures also show the accuracy of the pN approximations for the amplitudes 
of the three dominant modes. In each case, the amplitude of the {i, m) = (2, 2) mode 
is the most accurate, becoming worse for higher modes. This is to be expected, as the 
relative pN order to which the amplitudes are known decreases with increasing £. In 
particular, the (2, 2) mode is known to relative 3-pN order, while the (3, 2) and (4,4) 
modes arc only known to relative 2-pN order. Nonetheless, because of their far smaller 
magnitude, these higher-order modes actually have comparable absolute accuracy. 

The (3, 2) mode is particularly interesting. In the case, its amplitude is 

comparable to that of the (4,4) mode. However, in the S'^^ case, the (3,2) mode 
is far smaller until the merger. For most of the inspiral, the pN amplitude error is 
very large — being off by roughly 80%. Again, however, this error is relative; the pN 
approximation correctly predicts that the amplitude should be quite small in this 
case, because of a cancellation between the leading-order nonspinning and spinning 
components of the amplitude. 

Finally, we note that both waveforms can be hybridized to pN waveforms at 
frequencies of Mw sa 0.035, though these hybrids are not necessarily accurate enough 
to be useful in parameter estimation for detector-data analysis. As in reference [11], 
we can estimate the error in any hybrid by measuring the mismatch [71, 72, 73] 
between each pair of hybrids formed with the various approximants TaylorTl-T4; 
the error estimate is the maximum such mismatch. Using the Advanced LIGO high- 
power noise curve with no detuning [74] to do this measurement for the system, 
we find mismatches larger than 0.01 for total masses below roughly 40 M©. This 

the TaylorT2 approximant). That frequency never reaches the initial frequency of our simulation; 
the frequency plot "turns over" before that point. This is not unusual behavior for TaylorTS. For 
example, even in the nonspinning case, similar behavior is seen for mass ratio 10 : 1. Evidently, the 
series-inversion procedure is not particularly accurate. 
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means that, for any detected (SNR > 8) system with a lower mass, the uncertainty 
in these hybrids would be larger than the statistical uncertainty due to noise in the 
detector [75, 76]. For stronger signals or lower masses, more accurate pN waveforms 
and/or longer numerical simulations would be needed. For the S^'^ similar 
comparison would lead us to conclude that the hybrid is completely uncertain because 
of the bad behavior of the TaylorT3 approximant. If, on the other hand, we exclude 
that approximant as anomalously bad, we find that the hybrids are only accurate 
enough for parameter estimation above roughly 60 Mq. Still, it is possible that such 
hybrids would be accurate enough for detection purposes [77]. 

5. Conclusion 

The simulations discussed in this paper contain the most nearly extremal BBHs 
simulated to date. In our spin 0.97 simulation, we have found remarkably good 
agreement between the horizons' mass and spin evolutions and Alvi's analytic 
predictions, but we have found only moderately good agreement between the remnant 
hole's final spin and several analytic formulas, which suggests that these analytic 
formulas could be improved significantly using a set aligned and anti-aligncd nearly- 
extremal BBH simulations. 

We have found that the waveform from the 12.5-orbit anti-aligned case S^"^ 
exceeds the NINJA-2 accuracy requirements, while the waveform from the 25.5-orbit 
aligned case 5°+'' exceeds the NINJA-2 amplitude requirement but (because it is so 
long) fails to meet the NINJA-2 phase requirement (although it does meet the phase 
requirement easily if truncated to 5 orbits, the NINJA-2 length requirement). 

These results demonstrate the feasibility of applying waveforms from numerical 
simulations to gravitational-wave data analysis efforts when the holes have nearly 
extremal spins — a case previously inaccessible numerically but relevant astrophysically, 
given the evidence that nearly extremal black holes could exist. For example, 
waveforms such as those considered in this paper could be used in calibrating analytic 
template banks used for gravitational-wave detection searches. To pursue this goal, 
we plan to apply our methods for evolving nearly extremal BBHs to a large variety of 
BBH configurations, including unequal masses and spin precession. 

We have compared our numerical waveforms to several pN approximants, finding 
that the pN and numerical waveforms disagree at times well before merger. We 
also find that the pN approximants disagree with one another, indicating a large 
uncertainty in the pN approximations which leads to a large uncertainty in the 
resulting hybridized waveforms. Extracting the BBH parameters from detector data 
for systems with nearly extremal spins will require far longer numerical simulations, 
far more accurate pN waveforms, or a combination of the two. In the absence 
of improved pN waveforms, however, this implies that parameter estimation when 
the holes have nearly extremal spins could prove quite challenging, since the longer 
numerical simulations that would be necessary will come at high computational cost. 
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